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From  PS-splines  to  NURPS 


Paul  Dierckx  and  Joris  Windmolders 


Abstract.  A normalized  B-spline  representation  for  Powell-Sabin  (PS) 
spline  surfaces  is  extended  to  piecewise  rational  surfaces  (NURPS).  We 
investigate  the  adaptation  of  existing  algorithms  operating  on  B-splines 
to  this  more  general  case,  the  influence  of  weights  and  their  geometri- 
cal interpretation,  the  possibility  of  representing  planar  sections,  and  the 
conversion  from  rational  Bezier  to  NURPS  surfaces. 


§1.  Basic  Concepts 

1.1.  PS-splines 

Let  0 C R2  be  a simply  connected  subset  with  polygonal  boundary  Sft.  Let 
A be  a conforming  triangulation  of  Q having  n vertices  V)  with  coordinates 
( Ui,Vi ),  i = 1, . . . , n,  and  let  A*  be  a Powell-Sabin  (PS)  refinement  of  A (see, 
e.g.  [3]),  where  each  triangle  p € A is  divided  into  6 subtriangles.  A Powell- 
Sabin  (PS)  spline  is  a piecewise  quadratic  polynomial  with  Cl  continuity  on  Q. 
Dierckx  [1]  shows  how  to  calculate  a normalized  B-spline  basis  for  PS-splines: 

Definition  1.  A PS-spline  surface  has  a normalized  B-spline  representation 

n 3 

s(u,  = C:JS<  ( u > v ) € (!) 

i=l  j= 1 

where  c-,j  = ,(cfj , c?j , c?j)  are  the  B-spline  control  points  and  Bj(u,v)  are 
the  normalized  B-splines. 

This  representation  shares  a number  of  properties  with  tensor-product 
B-splines,  making  it  a powerful  tool  for  representing  surfaces  in  CAGD.  We 
summarize  the  most  important  properties  here.  For  details  we  refer  to  the 
original  paper  [1], 
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Fig.  1.  Domain  triangle. 

Property  1.  {B? (u,  v)}3^’2'^  is  a partition  of  unity: 
Bf(u,v)  > 0,  (u,  v)  G ft, 
£"=i£-=i^M  = i,  Mg  ft. 


Furthermore,  Bi}j(u,v)  is  nonzero  only  on  triangles  p G A having  V, \ as  a 
vertex: 


Property  2. 


= o, 


l 7^  i. 


(2) 


The  local  control,  affine  invariance  and  convex  hull  properties  follow  im- 
mediately. Linear  functions  can  be  represented  exactly.  In  particular,  we  will 
make  use  of  the  representations 

u = J2Yl  Ui’iBi  (u>  *')’  v = 

1 — 1 j — 1 t=l  j = 1 


Definition  2.  The  PS-triangles  U{Qit\,Q^2,  Qi,  3),  l = 1,. . . ,n  in  the  planar 


domain  have  as  vertices  the  B-spline  ordinates  Qi,j{Uij ,Vij),  j = 1,2,3. 

Consider  a domain  triangle  Pi,j,k{Vi,Vj,Vk)  G A with  its  PS-refinement 
(see  Figure  1).  Denote  the  Bezier  ordinates  as  svj,  v = i,j,k,l  = 1,2, 3, 4; 

(l,m)  G {(i,j),(j,k),(k,i)}  and  vi<jtk-  They  can  be  written  as 
unique  barycentric  combinations  of  the  B-spline  ordinates: 

Sy,t  — &v,l  Qv,  1 4"  Pv,l  Qv,2  4"  'tvj  Qv, 3i  (3) 

tl,m  = ^ l,m  $1, 2 4"  ^l,m  Sm,3i  (4) 

ul,m  = 3l,m  sl, 4 4"  m $m,4)  (5) 

C, ,j  ,k  — &i, 4 4"  Miq./c  Sjt4  4~  ^i,j,k  (fi) 


For  a given  PS-refinement  A*,  the  position  of  the  Bezier  ordinates  is 
fixed.  This  is  not  the  case  for  the  B-spline  ordinates.  The  following  lemma 
however  states  that  there  is  a restriction  on  the  B-spline  ordinates. 
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Lemma  1.  In  order  for  the  basis  functions  {Bj  (u,  i ;)}£z*’2’3n  to  constitute  a 
partition  of  unity  on  fl,  it  is  required  that  for  each  vertex  Vi,  i = 1, . . . ,n, 
the  PS-triangle  tj(Q»,i,  Qi,2,  Qi,3)  contains  the  Powell-Sabin  points,  i.e.,  the 
Bezier  ordinates  Sij,  l = 1,2, 3, 4,  of  any  domain  triangle  having  Vi  as  one  of 
its  vertices. 

There  is  a one-one  connection  between  the  barycentric  coordinates  of  the 
Powell-Sabin  points  at  vertex  Vi  with  respect  to  ti  and  the  value  of  the  basis 
functions  Bj(u,v),  j = 1,2,3,  and  of  their  derivatives  at  Vi,  e.g. 

Bj  (Vi)  = aM,  Bf  (Vi)  = A,i,  Bf  (Vi)  = 7i,i-  (7) 

Given  a PS-spline  surface  (1),  the  corresponding  Bezier  net  can  be  calculated 
efficiently  by  using  convex  barycentric  combinations  of  the  B-spline  control 
points  only: 

Property  3.  Applying  equations  (3)-(6)  where  the  ordinates  are  replaced  by 
control  points,  yields  the  corresponding  Bezier  net  of  the  surface. 

Finally,  via  the  concept  of  control  triangles,  the  B-spline  control  points 
give  us  valuable  insight  into  the  shape  of  the  surface: 

Definition  3.  The  control  triangles  are  defined  as  Tj(ci,i,  01,2,01,3). 

Property  4.  Each  control  triangle  T}(ci,i,ci, 3,01,3)  is  tangent  to  the  PS - 
surface  at  s (Vi). 


1.2.  NURPS 

The  Normalized  B-spline  theory  for  PS-surfaces  can  now  be  extended  to  a 
rational  scheme  just  like  tensor  product  B-splines  are  extended  to  NURBS. 
Referring  to  Figure  1,  we  use  the  boldface  notation  for  the  Bezier  points, 
e.g.  sv,j . Points  in  homogeneous  space  get  a h-superscript,  e.g.  s^j.  Their 

components  are  (s*f , s*f,  , s“,) . 

Definition  4.  A Non  Uniform  Rational  Powell-Sabin  (NURPS)  spline  surface 
has  the  form 


s(u,  v ) 


T,l=iH3j=iwi,jBJi(u,v) 


(u,v)  € fi, 


(8) 


where  cy  = (cy , c)V , c*j ) are  the  B-spline  control  points.  We  impose  that 
Wij  > 0 in  order  for  s(u,  v ) to  be  defined  anywhere  on  fl. 

If  u>ij  = 1,  i = 1, . . . , n,  j = 1,2,3,  then  (8)  reduces  to  (1).  The  following 
properties  are  readily  verified: 


48 


P.  Dierckx  and  J.  Windmolders 


Property  5. 


n 3 


s(u,v)  = #(«,«), 

1=1  j = 1 


(9) 


where 


and 


#(«>«) 


Wi'jBl{u,v ) 


SLi  E,=i 


(10) 


J 4>\{u,v)  > 0,  (u,i>)  € ft, 

1 E"=i E;3=i = i»  (w-^)  e n. 

Furthermore,  v)  is  nonzero  only  on  triangles  p g A having  Vi  as  a 

vertex. 


This  again  implies  the  local  control,  affine  invariance,  and  convex  hull 
properties. 

Property  6.  A NURPS  representation  (8)  is  the  3D-projection  in  Euclidean 
space  of  a 4D  PS-spline  in  homogeneous  space: 


s(u,v)  = (u’v)> 

*=1  j— 1 


(11) 

(12) 


§2.  Evaluation  and  Subdivision 

The  evaluation  of  s(u , v)  is  performed  in  two  steps: 

• First,  the  corresponding  rational  piecewise  Bezier  representation  is  cal- 
culated. 

• Then,  the  rational  de  Casteljau-algorithm  calculates  a point  on  this  ratio- 
nal piecewise  quadratic  Bezier  surface.  This  section  shows  how  to  perform 
the  first  step  in  a numerically  stable  way.  For  the  second  step,  we  refer 
to  Farin  [2],  Chapter  17.9. 


2.1.  In  homogeneous  space 

Formulae  (3)-(6)  can  be  applied  directly  in  homogeneous  space,  e.g. 


sv,l  = at>, I cv,l  +Pv,l  Cv,2  + 7v,i  cv,3 

- (<th’x  sh'y  th'z  \ 

~ >*V,1  >SV,l  J 

v!j,k  = Kj,k  Si)4  + Pi,j,k  Sy4  + Vi,j,k  Sk>4 

(htx  h.y  h.z  w \ 


(13) 


(14) 
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Projection  back  to  Euclidean  space  yields 


sv,i  = 


h.x 

SV,l 

„h>y 

Sv,l 

oW  '■ 

Sv,l 

’ 0*0 
sv,l 

Ch’z 

Sv,l 

Sv,l 


viJ,k  = 


h.x 

Vi,j,k 


n,h’V  „,h’z 


ViJ,k 


ui,j,k 


This  algorithm  has  a serious  drawback:  if  the  weights  vary  greatly  in 
magnitude,  the  coordinates  s£’(r,  v^Jk,  r = x,y,z  are  blown  away ; the  cal- 
culations don’t  operate  in  the  convex  hull  of  the  control  net  anymore,  and 
numerical  stability  is  endangered. 


2.2.  A rational  algorithm 

The  idea  behind  the  rational  de  Casteljau-algorithm  from  Farin  [2]  is  to  im- 
prove numerical  stability  by  rearranging  the  calculations,  avoiding  working  in 
homogeneous  space: 

sv,l  - av,l  wv,l  + Pv,l  Wv,2  + 7v,i  wv,3-  (15) 


Set 


- av,l  wv,l  - n 3 Pv  A Wv  2 n _ 7v,f  ®t>,3  v.  n 

«»,i  = ~w  P °>  A»,I  = —*3 > 0,  7v,(  = > o.  (16) 


°v,l 


*vA 


Then 

with 


®v,l  — ®v,l  Cv,X  T Pv,l  Cv,2  "l"  7v,i  CV(3 


&v,l  + Pv,l  4“  %,l  — 1. 


(17) 


(18) 


The  point  sv,i  is  a convex  barycentric  combination  of  cv,i,cVi2  and  cv,3,  so 
numerical  stability  is  guaranteed.  Likewise,  we  find 


tw  C W I 

l,m  ~ °ltm  % 2 ' €l,m  5m,3> 

uT,m  = &l,m  S™4  + €( ,m  sm,4> 

vf,j,k  = ^4  + SJ,4  + vi,j,k  *k,4> 


(19) 

(20) 
(21) 


C S™2  . _ el,rn  sm, 3 

Ol,m  — ,w  , eJ,m  — ,w  i 
t/l,m  * l}m 

c*  _ ^ l>m  Sf?4  _ el ,m  «m,4 

i €i,m  > 


^1,171 


\ ^ i,j,k  sT,4  - Pi,j,k  s^4 

*i,j,k  — i Pi>j,k  — > vi,j,k  — 


Vi,j,k 


Vi,j,k 


ui,j,k  $k,4 

dw  ’ 
Vi,j,k 


(22) 

(23) 

(24) 


^i,m  4”  6i,m  — 4"  ^ l,m  ~ ^ i,j,k  4”  Pi,j,k  4"  &i,j,k  — lj 


where 
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and  finally 


= $l,m  Si(2  “h  Q,m  Sm>3, 

(25) 

^ lfm  ®1,4  ^l,m  ®m,4  5 

(26) 

= ®i,4  4“  Sj}4  "b  V i,j,k  Sk,4* 

(27) 

All  formulae  are  convex  barycentric  combinations  operating  in  the  convex 
hull  of  the  B-spline  control  net.  After  having  computed  the  rational  Bezier 
representation,  Farin’s  rational  de  Casteljau  algorithm  can  be  used  to  evaluate 
the  surface  at  any  point  (u,v)  g fi. 


2.3.  Subdivision  on  uniform  triangulations 

The  evaluation  and  subdivision  of  spline  curves  and  surfaces  are  closely  related 
problems.  For  the  particular  case  of  a uniform  triangulation  A,  a subdivision 
scheme  for  PS-surfaces  has  been  derived  [4].  As  an  application,  it  was  shown 
how  a wireframe  of  the  surface  can  be  calculated  in  an  efficient  and  numeri- 
cally stable  way.  This  scheme  can  easily  be  extended  to  NURPS  on  uniform 
triangulations  again  using  Farin’s  technique  from  the  previous  section.  The 
details  are  omitted  here. 


§3.  Control  Planes 

Recall  that  the  NURPS  representation  inherits  the  convex  hull,  affine  invari- 
ance, and  local  control  property  from  the  normalized  B-spline  representation. 
This  section  adds  the  tangent  property  to  the  inheritance  list,  and  shows  how 
the  rational  representation  allows  for  more  flexibilty  when  designing  surfaces. 


3.1.  Tangent  property 


Referring  to  the  locality  of  the  B-splines  (2),  it  is  easy  to  verify  that  the 
evaluation  of  s(u,  v ) and  its  derivatives  at  vertex  Vi  yields 


s{Vi ) = ait  1 Ci,i  C;,2  + 7i,l  cii3, 


ds(V j) 
du 

ds(V j) 
dv 


Ci,i  + et-,2  c;,2  + e*,3  ci>3, 
Ci,l  + dj, 2 <7,2  + d^3  cii3, 


(28) 

(29) 

(30) 


for  some 

ei,l  + ei,2  + ei,  3 = diti  + dit  2 + d,-,3  = 0. 

It  follows  that  the  control  triangle  at  V)  is  tangent  to  the  surface  at  s(Vj),  i.e. , 
any  point  p in  the  tangent  plane  is  a barycentric  combination  of  the  control 
points  Ci,i,Cii2,cii3: 


P = s(Vi)  + a 


dsjVj) 

du 


+ b 


ds(V j) 
dv  ’ 


a,  b G IR. 


This  is  illustrated  in  Figure  2 (left). 
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Fig.  2.  NURPS  surface  and  its  control  planes;  local  planar  effects. 


3.2.  Shape  parameters 


Farin  [2]  introduces  the  concept  of  shape  parameters  with  respect  to  rational 
Bezier  curves.  A geometric  handle  allows  the  designer  to  influence  the  shape 
of  the  curve  in  a predictable  way,  rather  than  requiring  the  input  of  numbers 
for  the  weights.  In  the  same  work,  it  is  stated  that  this  property  does  not 
carry  over  to  rational  Bezier  surfaces  on  triangles,  but  shape  parameters  can 
be  defined  for  NURPS._ 

Recall  that  (aV]i, /3v,i>7t),i)  are  the  barycentric  coordinates  of  sv>i  with 
respect  to  control  triangle  rv(cVli,cV)2,cV)3).  From  (16)  it  follows  that  sv,x 
can  be  moved  within  Tv  to  a new  location  (&v,i,pv>i,%ti),  while  keeping  its 
weight  s™,  constant.  The  corresponding  PS-weights  are  found  immediately 
as 


Wv,l  = 


WV}2  = 


0V,1 


Wv,  3 = 


%,1<1 

7v,i 


(31) 


This  shows  how  (a^i, j3v,i,  7w,i)  can  be  used  as  shape  parameters. 


3.3.  Planar  sections 


Definition  5.  Let  [tj, t2,  • • • ,t„]  denote  the  convex  hull  of  the  3D  points 
tl, t2, ■ • • , tn. 

Definition  6.  Let  5(A)  denote  the  image  of  a subset  A C under  (8). 

Definition  7.  Let  r(a,b,  c)  denote  the  Bezier  subtriangle  in  the  domain  plane 
with  vertices  a,  b and  c. 

If  the  control  triangles  of  adjacent  vertices  Vi,Vj,Vk  are  chosen  to  be 
coplanar,  then  the  surface  section  5 {pi,j,k  {Vi,  Vj,  V*,))  will  be  in  the  same 
plane,  as  a consequence  of  the  convex  hull  property.  However,  using  the 
weights  in  the  NURPS  representation,  it  is  possible  to  achieve  more  local 
planar  effects. 

The  rational  evaluation  algorithm  from  Section  2.2  reveals  that  for  = 
Wj, 2 = Wifi  = w — > oo,  i € (1, . . . , n)  and  referring  to  Figure  1,  the  following 
holds  on  the  domain  triangle  Pi,j,k(Vi,Vj,Vk)'. 


*ij  = 


fii,j 


Sij  + - 


‘■■’7,3 


-Si  2 + 


C*,J  + w 


-Sj,3- 


w 


(32) 
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Fig.  3.  Bezier  triangle. 


Thus, 


lim  tjj  = Sj  2 ■ 

to— *00 


Likewise,  for  the  other  Bezier  points  of  r we  find 

Si.i  = aiti  Ci,i  + pi, i cii2  + 7 i,i  <7,3,  / = 1, 2, 4, 
lim  u;j  = lim  viJik  = sii4. 

to— too  to— +oo 


Consequently, 

S (t  (Sj,1,  t{,j , Vi  J,k))  — [®i,  1 , ®i, 2)  ®i,4 ] • 

Similar  reasoning  on  the  other  Bezier  subtriangles  shows  that 

S(T(Si,i,Vij,k,tk>i))  = [Sj,i , Sj^,  Sjt3]  , 

S {T  {ti,j  > si,l)  — [Si,2,Sj,4]  , 

S (r  {tk,iiVi,j,k , Sk  i ))  = [s j,3 , Si,4]  , 

S {t  (,Vi,j,k,  Sj,\,tj,k))  — [Si,4]  , 

S (t  (,Vi,j,k,tj,kt  sk, l))  — [si,4]  i 

and  therefore, 


(33) 

(34) 

(35) 


S(piJ,k(Vi,Vj,Vk))  = [Si,l » Sj,2 , S|,4 , S|,3]  C [cj,x , Cj,2 , Ci,3] 

The  latter  image  is  a planar  surface  section.  Figure  2 (right)  shows  some 
NURJPS  surface  with  very  large  weights  at  a vertex. 


§4.  Conversion  from  Rational  Bezier  to  NURPS  Representation 

Suppose  we  are  given  a rational  quadratic  Bezier  surface  on  one  domain  tri- 
angle (see  Figure  3)  p(si,i,  sj}i,sk,i) 

b{u,  v)  = Y uhits),  (36) 

*l+*2+*3  = 2 

where  *1,12,13  > 0,  ( u,v ) € p and  (<i,t2>^3)  are  the  barycentric  coordinates  of 
( u , v)  with  respect  to  p.  In  this  section  it  is  shown  how  a NURPS  representa- 
tion 

S(U,V)=  Y Y Cl,mBin(U'V)  (37) 

l~i,j,k  m=l 

of  the  given  surface,  for  a specific  choice  of  the  PS-triangles,  is  immediately 
obtained.  To  simplify  the  notation,  the  surfaces  are  considered  in  homoge- 
neous space. 
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Fig.  4.  Subdivision  at  and  t^j. 

Lemma  2.  Suppose  we  are  given  a triangle  t (Vi,  Vj,  Vjt)  with  barycenter  z. 
IfWi  denotes  the  midpoint  of  the  side  opposite  to  VJ,  then  (z  + Vi)  / 2 is  the 
barycenter  of  the  triangle  t (Vi,  Wm,  Wn)  ,l,m,n  6 {i,j,  k]  ,1  ^ rn  ^ n. 

The  construction  of  the  NURPS  representation  relies  on  the  de  Castel- 
jau-algorithm  for  Bezier  triangles  (see,  e.g.,  [2]).  Subdivision  at  the  barycenter 
of  p and  at  the  midpoint  of  edge  Sj.iSj.i  (see  Figure  4)  yields  the  new  Bezier 
points 


bi!o,o  = g (^2,0,0  + ki.i.o  + bj  01) , 

(38) 

boii.o  = 3 (b  1,1,0  + bo)2,o  + bS,i,i)  , 

(39) 

^oio.i  = 3 (^i,o,i  + ^0,1,1  + bo, 0,2)  > 

(40) 

bo, 0,0  = 3 (^i!o,o  + ^o!i,o  + bo!o,i)  > 

(41) 

^1(0,0  = 2^2,o,o  + bf,i,o)) 

(42) 

do’l.0  = 2 (^1.1,0  + ^o,2,o)> 

(43) 

jl,/i  /i  lth  1 \ 

ao,o,i  — 2^di.o,o  Doti,oJi 

(44) 

j2,Jl  1 / jl,/l  _L  \ 

ao,o,o  ~~  2 ^aM>o  ^ ao,i,oJ* 

(45) 

After  subdivision  of  the  two  remaining  edges,  the  6 subtriangles  thus  obtained 
constitute  a PS-refinement  of  p,  say,  with  interior  point  Vij^  and  edge  points 
ti,j>tj,k>tk,i  (see  Figure  5,  left).  A NURPS  representation  of  the  given  surface 
on  this  PS-refinement  is  easily  obtained  (see  Figure  5,  right).  Set 

Qv,  1 ^ Jj 

Qif 2 — Qj',3  = ti,j  j Qj,2  = Qfc,3  = tj,ki  Qfc,2  — Qi,; 

3 — 

„h  _ lJi  f.h  _ l/i  h _ 

ci,l  ~ °2,0,0’  ci,2  — D1,1,0>  ci,3  — D1,0,1> 

rh  _ uh  „h  _ lIi  h _ 

cj',l  ~~  °0,2,0’  cj, 2 — D0,1,1>  cj, 3 — D1,1,0> 

_ u.h  „h  _ uk  nh  _ uft 

cfc,l  ~ D0,0,2>  cfc, 2 — D1,0,1>  ck, 3 — D0,1,1‘ 


and 
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Fig.  5.  PS-refinement. 


Then  by  Lemma  2,  it  follows  that  the  PS-points  svj,  l = 1, 2, 3, 4 are  inside  the 
PS-triangle  tv,  v = i,j,k.  Now  recall  formula  (13)  for  v = i,j,k,  and  l = 4, 
and  formula  (14),  with  aVi 4 = /3„,4  = 7^,4  = 5,  resp.  A = fiij.k  = ux,j,k  - 
| , in  order  to  calculate  the  corresponding  Bezier  points  of  this  NURPS  surface. 
It  turns  out  that  these  equations  are  exactly  the  same  as  the  subdivision 
formulae  (38)-(41).  Likewise,  since  in  (3)-(5) 

(®v,2i  fiv,2x  7u,2)  ( 2 ) 2 ’ ^)’ 

1 1 

(^v,3i  'Tv, 3)  (^  > 2)' 

1 1 

= (^> 

for  v = i,j,k  and  ( l,m ) € {(w),C?,  A:),(fc,i)}  , similar  reasoning  shows  that 
calculating  the  corresponding  Bezier  net  of  (37)  exactly  yields  the  Bezier  net 
of  (36)  after  the  proposed  subdivisions.  Hence,  b(u,v)  = s(u,v)  on  p. 
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